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ABSTRACT 

O , The dynamics of a viscous accretion disc subject to a slowly varying warp of large 

' amplitude is considered. Attention is restricted to discs in which self-gravitation is 

(— i , negligible, and to the generic case in which the resonant wave propagation found in 

' inviscid Keplerian discs does not occur. The equations of fluid dynamics are derived 

in a coordinate system that follows the principal warping motion of the disc. They are 
reduced using asymptotic methods for thin discs, and solved to extract the equation 
governing the warp. In general, this is a wave equation of parabolic type with non-linear 
dispersion and diffusion, which describes fully non-linear bending waves. This method 
, generalizes the linear theory of Papaloizou & Pringle (1983) to allow for an arbitrary 

£NJ ' rotation law, and extends it into the non-linear domain, where it connects with a 

generalized version of the theory of Pringle (1992). The astrophysical implications of 
this analysis are discussed briefly. 
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1 INTRODUCTION 



OO 
G\ 

Pl,: 

6 

■ There are many situations in astrophysics in which accretion discs are believed to be non-planar, either because their profile is 
| observed directly (e.g. in the nucleus of the galaxy NGC 4258; Miyoshi et al. 1995) or in order to explain phenomena such as 

• ^ ■ the periodic modulation of light curves (e.g. in the X-ray binary Her X-l) or the precession of jets (e.g. in the X-ray binary SS 
rS ' 433). In many - though not all - applications the effects of self- gravitation can be neglected and the evolution is dominated by 

■ viscous fluid dynamics. The earliest attempts to derive equations governing the time-dependence of a slowly varying warp of 
small amplitude (Petterson 1978; Hatchett, Begelman & Sarazin 1981) suggested a simple diffusive evolution but were shown 
by Papaloizou & Pringle (1983) to be incorrect at a fundamental level, being inconsistent with the conservation of angular 
momentum. 

Papaloizou & Pringle (1983) demonstrated that the problem is complicated because of the importance of a resonance in 
thin discs that are both Keplerian (or very nearly so) and inviscid (or very nearly so). It is therefore necessary to distinguish 
between the generic (non-resonant) case and the resonant case which occurs when both 



n- 



<, H/r and a <, H/r. (1) 



Here Q is the angular velocity, k is the epicyclic frequency, H/r is the ratio of the semi-thickness of the disc to the radius, 
and a is the dimensionless viscosity parameter (Shakura & Sunyaev 1973). Unlike most resonances in discs, which occur at 
specific radii, this one, which is effectively the m — 1 inner Lindblad resonance, is likely to occur either everywhere in the disc 
or nowhere. Current estimates suggest that the resonant case is not relevant to discs in X-ray binaries or active galactic nuclei, 
but may apply in protostellar discs. However, the resonance is delicate and could be destroyed by the effects of self-gravitation, 
magnetic fields, or turbulence. 

Papaloizou & Pringle (1983) derived the first consistent equation governing a warp of small amplitude. They considered 
the (non-resonant) case of a Keplerian disc with a significant viscosity, treating the warp as a slowly varying disturbance with 
azimuthal wave number m = 1, using linear Eulerian perturbation theory. The resulting equation for the warp is a complex 
linear diffusion equation, which has been used in applications (e.g. Kumar & Pringle 1985). 
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Papaloizou & Lin (1994)0 considered the case of an inviscid disc. Again, linear Eulerian perturbation theory was used, 
and the authors assumed the warp to be a slowly varying normal mode of the disc. In the non-Keplerian (non-resonant) case 
the warp obeys a dispersive linear wave equation, while in the Keplerian (resonant) case the warp obeys a non-dispersive linear 
wave equation. By considering the effect of a small viscosity on the inviscid modes, Papaloizou & Lin (1994) connected this 
theory with that of Papaloizou & Pringle (1983), showing how the transition occurs between wave-like and diffusive behaviour 
in Keplerian discs when a ~ H/r. This theory has also been used in applications (e.g. Papaloizou & Terquem 1995). 

The major uncertainty in these theories is whether the linear analysis is valid for warps of a sufficient amplitude to be 
observable. The linear Eulerian perturbation theory is formally valid only when the tilt angle f3(r, t) of the warp satisfies \/3\ <C 
H/r. However, because of the special nature of the mode involved, which differs little locally from a rigid tilt /3 = constant, 
one suspects that the linear theory might be valid for larger warps. In fact, a more appropriate measure of the amplitude of 
the warp is \d/3/d\nr\. The Eulerian perturbation theory offers little information about possible non-linear effects, although it 
does predict that, in the resonant case, an amplitude \d(3/d\nr\ ~ H/r would result in horizontal shearing motions in the disc 
comparable to the sound speed, which are expected to be unstable (Kumar & Coleman 1993; Gammie, Goodman & Ogilvie, 
in preparation). 

Meanwhile, Pringle (1992) developed a different approach in which the forms of the equations governing a warped viscous 
disc are derived simply by requiring mass and angular momentum to be conserved, but without reference to the detailed 
internal fluid dynamics of the disc. In this scheme, neighbouring rings in the disc exchange angular momentum by means of 
viscous torques which are of two kinds. One kind of torque (associated with a kinematic viscosity coefficient v\) acts on the 
differential rotation in the plane of the disc, and leads to accretion; the other kind (associated with 1/2) acts to flatten the 
disc. The equations derived are 

^ + ~(rEe r ) = (2) 
at r or 

for the surface density E(r, t), and 

* (E^IM) + I » (Eft.r'M) = I * U^fi) + 1 » (h^r^f) (3) 
at r or r or \ dr J r or or J 

for the angular momentum, in the absence of external torques. Here v r (r,t) is the mean radial velocity and £(r,t) is the tilt 
vector, which is a unit vector parallel to the orbital angular momentum of the ring. From this follow the equations 

^ d(r 2 D.) 13 / r s dQ\ , ^ . d£ 2 

dr r or \ dr J or 

for the component of angular momentum parallel to £, and 



(4) 



d£ ( dlnfi\ d£l 13 /, ^ 3 „d£\ , ^ 2r, 

at \ dr J or J ror V' or J * 



(5) 



for the tilt vector. 

This approach appears to be valid for warps of large amplitude and would therefore represent an advance on the linear 
theory. It has been used in several applications, in particular to identify the radiation-driven instability (Pringle 1996) and 
to explore its linear theory (e.g. Maloney, Begelman & Pringle 1996) and non- linear evolution (e.g. Pringle 1997). However, 
because the equations of Pringle (1992) are derived somewhat heuristically, without reference to the detailed internal fluid 
dynamics of the disc, some doubts remain over the validity of this method. The obvious questions to be addressed are whether 
any internal degree of freedom of the rings has been neglected, whether the interaction between neighbouring rings is purely 
of the assumed form of viscous torques, and whether there are any non-linear fluid-dynamical effects that might limit the 
amplitude of the warp. It is therefore important to investigate whether the equations of Pringle (1992) can be derived ab 
initio from the three-dimensional fluid-dynamical equations, and to understand how they connect to the previously established 
linear theory. These are the aims of this paper. 

The analysis is organized as follows. I first define a coordinate system that follows the principal warping motion of the 
disc, and derive the basic results necessary for vector calculus (Section 2). I then present the exact forms of the equations of 
fluid dynamics in this coordinate system (Section 3). The equations are reduced using an asymptotic analysis for thin discs 
(Section 4). I then apply separation of variables to solve the fully non- linear problem except for the determination of three 
dimensionless coefficients from the solution of a set of ordinary differential equations (ODEs) (Section 5). These are evaluated 
in Section 6. Readers interested only in the interpretation and application of this work may proceed to Section 7, where the 
principal results are summarized and discussed further. 



* See also Papaloizou & Lin (1995), in which self-gravitation is included. 
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/(r,<)=(sinjScos7,sin(Ssin7,cos|8) 




Figure 1. Warped spherical polar coordinates (r,0,<j>). The axis of each sphere r = constant is tilted to point in the direction of the 
unit vector £(r, t). 

2 DEFINITION OF WARPED COORDINATES 

The non-linear fluid dynamics of a warped disc is most naturally described in the case of a spherically symmetric external 
potential. Then a flat disc has no preferred plane of orientation, and must possess a zero-frequency mode consisting of a 
rigid tilt of any amplitude. Continuous with this mode, there exist bending waves, with azimuthal wave number m = 1 in 
linear theory, which vary on a time-scale long compared to the orbital time-scale, and on a length-scale long compared to the 
thickness of the disc. This is in contrast to the other modes of the disc (Lubow & Pringle 1993) and motivates the following 
analysis. 



2.1 Coordinates and basis vectors 



Let (x,y,z) be Cartesian coordinates in an inertial frame of reference, and define the spherical radial coordinate r = (x 2 + 
y 2 + z 2 ) 1 / 2 . On each sphere r = constant, define the usual angular coordinates (6, (j>), but with respect to an axis that is tilted 
to point in the direction of the unit vector £(r,t) (Fig. 1). It is intended that the disc matter on each sphere will lie close to 
8 — tt/2. The tilt vector can be described by Euler angles /3(r,t) and 7(r,t): 



i = sin P cos 7 e x + sin /3 sin 7 e y + cos fi e z . 

In detail, warped spherical polar coordinates (r, 0, 



are defined by 



x 

y 

z 

where 



M(r,M,*) 



M(r,M,t) = 



cos 7 
sin 7 




— sin 7 
cos 7 




cos p 


. — sin P 



sin P 


cos p. 



cost 
sin (, 




— sin i 

cos <j) 




sin t 


cos ( 



COS U 



— sin 6 



(6) 



(7) 



(8) 



is a composition of four orthogonal linear transformations. For the functions P(r,t) and 7(r, t), a prime and a dot will denote 
differentiation with respect to r and t, respectively. 

Although (r, 6, (f>) are not orthogonal coordinates, it is appropriate to refer vectors to a natural orthonormal basis 
{e r , eg, e^} defined by 
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" e r " 


e y 


= M 


eg 






. &4, . 



The vector e r is normal to the sphere r = constant, while the vectors eg and e^ are tangent to it. The components of an 
arbitrary vector F then transform according to 



-F x - 




-F r - 


Fy 


= M 


F e 


.F„. 




If J 



2.2 Calculus in warped coordinates 

The matrix M is orthogonal. Its decomposition in equation (^) makes it easy to compute the differential 
dM = MA, 

where A is an antisymmetric matrix, with components given by 
A12 = — d9 — cos 4>&f3 — sin j3 sin (j>drf, 

An — — sin 9 d<j) + cos 9 sin cf> d/3 — (cos /3 sin 9 + sin /3 cos 9 cos cf>) d7, 
A23 = — cos 9 d<f> — sin 9 sin d/3 — (cos/3 cos — sin /3 sin cos 0) d7- 

The coordinate differentials are therefore related by 

dr 



"da;" 




dy 


= M 


.d0. 





r/3' cos + rj 1 sin /3 sin 



rd9 + r cos d/3 + r sin /3 sin </> d7 
r sin 9A(f> — r cos (9 sin </> d/3 + r (cos /3 sin + sin /3 cos 6 1 cos <f>) d'y 

so the Jacobian matrix of the transformation is 
'dx/dr dx/d9 dx/dcj>' 
dy/dr dy/d8 dy/dcj> = M 
.dz/dr dz/ 89 dz / d<f> . 

with determinant 
J — r 2 sin 9, 

as for ordinary spherical polar coordinates. 
The gradient operator takes the form 



. —r/3' cos 9 sin </> + vy' (cos /3 sin 9 + sin /3 cos cos <j>) r sin 9 . 



e r T> + ee -de + e 4 



r sin 9 



where 



V 



d r — (/3' cos + 7' sin /3 sin </>)<9e — [— /3' cos sin <j> + 7' (cos /3 sin + sin /3 cos cos c 



sin# 



(),<, 



is a modified radial derivative. 
Finally, the condition 

~de x ' 
de v = 
_de z . 

yields the differentials of the unit vectors in the form 

de r = eg (d9 + cos </> d/3 + sin /3 sin d7) + e^ [sin # d</> — cos sin </> d/3 + (cos /3 sin 9 + sin /3 cos 9 cos </>) d7] , 
deg = — e r (d# + cos cj> d/3 + sin /3 sin </> d7) + e$ [cos # d</> + sin 9 sin d/3 + (cos /3 cos # — sin /3 sin # cos <j>) d'y] . 
de^ = — e r [sin d</> — cos 9 sin </>d/3 + (cos /3 sin # + sin /3 cos 9 cos </>) d7] 

— ee [cos 6> d</> + sin 9 sin </> d/3 + (cos /3 cos # — sin f3 sin # cos <f>) d'y] . 

In particular, 



Ve r 
Ve« 
Ve^, 



1 1 

- egeg H — erf,e<A, 
r r 

1 cot / 
eee,- H e^e^ H : — - e r e^, 



1 cote* 

e^e r e^eg 



r sin 9 
f 

1 

r sin 



(ii) 

(12) 
(13) 
(14) 

(15) 

(16) 
(17) 

(18) 
(19) 

(20) 

(21) 
(22) 

(23) 

(24) 
(25) 
(26) 
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where 

/ = r/3' sin <j> — rj sin /3 cos <f). (27) 

The advantages of this method over that of Petterson (1978) are to be emphasized. Petterson (1978) introduced a warped 
coordinate system based on tilted cylinders, rather than spheres. This is less satisfactory because the resulting coordinate 
system has a non-trivial Jacobian determinant and, in fact, breaks down for large-amplitude warps. The present method does 
not suffer from these defects. Moreover, it can be understood using elementary vector calculus rather than requiring advanced 
differential geometry. 



2.3 Absolute and relative velocities 

Consider a fluid element with coordinates x(t), y{i) and z(t). The velocity components are 



~u x ~ 




'dx/dt- 




' u r ' 


Uy 




dy/dt 


= M 


Ug 


-U z _ 




.dz/dt. 







where 







Ug 









dr/dt 



r(&0/&b) + r(d(3/dt) cos (j> + r(dj/dt) sin /3 sin <f> 
r sin 0(dcj>/dt) — r(df3/di) cos 6 sin (j> + r(d7/d£)(cos /3 sin 6 + sin f3 cos 9 cos ( 



Here, d/dt stands for the Lagrangian time derivative 
D = (d t ) x ,y,z + u x d x + u y dy + u z d z . 



(28) 



(29) 



(30) 



The components (u r ,ug,u<p) are the absolute velocity components, in the sense that u — u r e r + ug eg + u^ is the velocity 
vector as measured in the inertial frame. It is intended, however, that the warped coordinate system itself will follow the 
principal motion of the fluid, other than its orbital angular velocity. The additional motion with respect to the warped 
coordinate system is described by the relative velocity components 



(31) 







dr/dt 


ve 




r(d0/dt) 


-V4,- 




_rsm6(d<j)/dt) 



The Lagrangian time derivative takes the form 

D = (dt)r,e,A + v r d r + -de A "^r—^ds. 

r rsmO 

The absolute velocity may then be written 
u = u r e r + ue eg + u$ e$, 
with 



Ur 
Ug 

U<t> 



vg + r(D/3) cos <f) + r(Dy) sin (3 sin <f), 

v$ — r(D/3) cos 6 sin <f) + r(Dy)(cos (3 sin9 + sin (3 cos 6 cos ( 



(32) 
(33) 

(34) 
(35) 
(36) 



3 FLUID DYNAMICS IN WARPED COORDINATES 
3.1 Basic equations 

The equations governing the dynamics of a compressible, non-self-gravitating fluid are the equation of mass conservation, 
Dp = -pV-u, (37) 
the adiabatic condition, 

Dp=-r P V-it, (38) 
and the equation of motion, 

pQu = -Vp - pV$ + V- [/jVu + /i(Vit) T ] + V [(/i b - §A0V-it] + F. (39) 
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Here p is the density, p is the pressure, V is the adiabatic exponent, $ is the external gravitational potential, p and p b are 
the shear and bulk viscosities, and F is an arbitrary external force. 

In using these equations, it is assumed that the agent responsible for angular momentum transport in accretion discs can 
be treated as an isotropic viscosity in the sense of the Navier-Stokes equation. In the case of magnetohydrodynamic turbulence, 
it is not known how accurate this assumption is, nor whether any alternative constitutive relation for the turbulent stress 
would be better. Although bulk viscosity is not usually discussed in the context of accretion discs, it is a standard component 
in the compressible Navier-Stokes equation and is included here for completeness. It is, however, assumed that the viscosities 
are only dynamically, rather than thermodynamically, important. A more sophisticated treatment would include the thermal 
effects of viscous dissipation and radiative transport. 

Using the information given in Section 2, it is possible to derive the form of these equations in warped coordinates. Either 
absolute or relative velocity components are used below, whichever gives the most compact form of the equations. First, note 
that V-u may be written in two equivalent ways: 



1 2 1 1 

V u = —T>(r u r ) H : — -de(u e sm9) H : — ^d^us 

r 2 r sin a r sin 9 

1 2 1 1 

= —zdrir v r ) H : — -dg(ve sm8) H : — tO^vj,. (40) 

r 2 rsmO r sin 

The familiar form of the latter expression is a consequence of the fact that the Jacobian determinant of the warped coordinate 
system is equal to that of ordinary spherical polar coordinates. This deals with equations (^) and (^). Equation ( |39| requires 
a lengthy calculation, especially to evaluate the viscous terms, although much of their complexity is already present in ordinary 
spherical polar coordinates. The result is 

2 2 N 

Du r - Hi - ^± ) = -Dp - pD<$> + F r 
r r 



+V \(p b + i^)V-ul + \v(pr 2 Vu r ) + 1 d e (psm9d e u r ) + 1 d^pd^Ur) - ^-T>(p,r) 
" v rsinfi r 2 sm 9 r J 

-^d 6 (u e sine) + (dep)V ( - ^fLfct* + M V + f ( J*^ - U ,Q 6 ,\ , (41) 

rsmS \ r J r" sm 9 sin # \ r J r 2 sin 9 \ sm 9 J 

p I Due + UrUe HL_ [uj, cos 9 + rCDS) sinr/> — r(Dy) sin 8 cos 6] X = ~-d g p — -9 e $ + F g 

I r r sin 9 J r r 

+ -d e I (jib + hfi)V-u] + \v{pr 2 Vu e ) + 1 d e (psm9d e ue) + 1 , d^pd^ug) 
v v r 2 sm 9 r 2 sm 9 

s u e co%9 . T>{pr 2 ) 9gp dg(psin 2 9) d^p ( \ 

^V(pr) . dg(psm9) H y —^-d e u r Vu r 9 ■ 3 n d 4> u * + 

r 2 r 2 sin 9 r d r r 2 sm 9 r 2 \ sin 9 ) 

1 -n ( ^ r l u 4> \ -r,„. vf u e W)\ 

r 2 \ sin^ J r sine r^sin 2 ^ ( ' 

p(dm + VzO± -i HL^ [u^cosfl + r(D/3)sin0 - r(D7) sin/? cos 0]) = 5— r—^d^ + Fa 

I r rsm.9 J rsinfl rsmff 

H r— [(Mb + |m)V-m] + -^©(^ir 2 ©^) + „ 1 dg(psin9 dgu^) + 1 „ d^pd^u^) 

r sin 9 v r 2 sin # sm & 

u <*-rv ^ ^cos6» . V{pr 2 ) d^p d B {psm 2 9) d^p / 119 \ 

V sm 6* / r sin 9 r 2 sm S 

These equations must be solved in conjunction with equations (pS4|)-(|36|) defining the absolute velocity components and 
equation (132) defining the Lagrangian time derivative. 



3.2 Generalized form of the angular momentum equation 

The equations above are capable of describing any flow, with /3(r, t) and 7(7", i) being arbitrarily prescribed (differentiable) 
functions. In order to derive meaningful dynamical equations for /3 and 7, one must impose suitable constraints on the flow. 
For a thin disc, the fluid must remain close to the surface 9 = n/2, and this condition is enforced by the asymptotic analysis 
presented in Section 4. It will be convenient to work with the dimensionless complex variable 

ip = t[>r + = r(/3' + ij'sin8). (44) 
Then \ip\ is a measure of the amplitude of the warp, with 
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= r 



d£ 

dr 



(45) 



It will be found that, although equation (g) for the surface density can easily be derived from the equations of this section 
in the limit of a thin disc, it is impossible to obtain equations (^) and Q in a similar manner. One must allow for a more 
general form of the torque between neighbouring rings, and attempt to derive an angular momentum equation of the form 



(46) 



|;(Er 2 ^) + ±|-(£*vr 3 fi£) = ~ (QiTr 2 Q 2 £) + l±( Q2 I r ^ 2 ^-) + 1°. ( Qi Tr^l x 
at r or r or y ' r or \ or J r Or \ Or 

Here I(r, t) is the azimuthally averaged second vertical moment of the density, which in cylindrical polar coordinates would 
be written 

<-2ir 



1 

2~7 



2d(j>, 



where T - 



pz dz, 



(47) 



and is an important dynamical quantity in the theory of bending waves. The coefficients Qi, Q 2 and Q3 are dimensionless 
quantities to be determined, which will be found to depend on the rotation law and the shear viscosity, and also, in the 
non-linear theory, on the adiabatic exponent, the bulk viscosity and the amplitude of the warp. From this follow the equations 



Sv r (r 2 n)' 



ld_ 

r dr 



Q 2 lr 2 Q 2 



dt 2 
dr 



for the component of angular momentum parallel to £, and 



> 2 fi ( 



d£ _ m 

dt r dr 



„ „ 2 d£ Id 
= QilrQ 2 — + - — 
dr r dr 



dr 



Q 2 Ir 2 Q. 2 



dl 2 „ Id 
or r dr 



(Q 3 lr 3 n 2 £ x g) 



for the tilt vector. The correspondence with /3 and 7 is made through the relation 
t = - e '\e=«/2> 

which can be differentiated using equation (j^f. Thus 
■£v r (r 2 n)' 



^d r (Qiir 2 Q 2 ) - Q 2 m 2 \i>\ 2 , 



Er 2 f2(/3 + Vr/3') 



X 



= QrJrfi 2 /3' + ~dr [Tr 3 ^ 2 {Q 2 - Q37' sin/?)] - lr 2 Q 2 -/' cos /3(Q 27 ' sin/? + Q 3 /3'), 

and 

Er 2 fi(7 + «r7')sin/? = Y 

= QiJrn 2 7'sin/?+ -d r [lr 3 n 2 (Q 2l ' sin P + Q 3 P')] + 2> 2 fi 2 7 ' cos P(Q 2 P' - Qal' sin 0). 
It will be useful to consider the complex combination 

X + IY = Q 1 lQ. 2 ii + -(d r + W cos/3)(Q 4 2:r 2 fiV), 
r 

where Q 4 = Q 2 + iQz- 



(48) 

(49) 
(50) 
(51) 

(52) 

(53) 
(54) 



4 NON-LINEAR BENDING WAVES IN A THIN DISC 
4.1 Asymptotic expansions 

Consider a thin disc in a spherically symmetric potential <&(r), subject to a slowly varying, unforced warp of large amplitude. 
Let the small parameter e be a characteristic value of the local angular semi-thickness of the disc. Then, to resolve the internal 
structure of the disc, introduce the stretched vertical coordinate £ by 

= ~ - <• (55) 

so that £ = O(l) within the disc. All quantities vary not on the fast, orbital time-scale, but on the slow time-scale characteristic 
of both viscous and dynamical evolution of the warp. This is captured by the slow time coordinate 

T = e 2 t. (56) 
[It is implicit that units are chosen such that the radius of the disc and the orbital time-scale are O(l).] Then set 

(3(r,t) 1 * /3(r,T), (57) 
7 (r,t) 1 * 7(r,T). (58) 
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A dot will now denote differentiation with respect to T. For the density and pressure, introduce the scalings 

p(r,$,4>,t) = e 3 [po(r,0,C,r) + e Pl (r,0,C,r) + O(e 2 )] , (59) 

p(r,M,t) = e 3+2 [po(r,0,(,T) + c Pl (r,0,(,T) + O(e 2 )] , (60) 

where s is a parameter which should be positive if the self-gravitation of the disc is to be negligible, but otherwise has no 
effect on the following analysis (cf. Ogilvie 1997). For the velocities, set 



v r (r,6,<t>,i) 



ev rl (r, </>, C, T) + eV 2 (r, 0, <, T) + 0(e 3 ), 



v e {r,e,<t>,t) = ev ei (r,<p,(,T) + e 2 V02(r,<l>X,T) + O(e 3 ), 



(61) 
(62) 
(63) 



iv(r,0,0,t) = rfl(r) sin 61 + e^i(r,0,C,T) + e 2 ^ 2 (r, </>, C, T) + C>(e 3 ). 

Finally, for the viscosities, assume 

fl(r,9^,i) = e s+2 [Mo(r,0,C,T) + e Ml (r,^,C,T) + O( e 2 )] , (64) 

p b (r,6\</>,t) = e s+2 [Mbo(r,0,C,r) + ^bi(r,0,C,r) + O(e 2 )] . (65) 

This is the correct scaling for an a- viscosity (Shakura & Sunyaev 1973) with a = 0(1), which includes the possibility of small 
or zero a unless resonance occurs. Note that, for this non-linear warp, all quantities other than the orbital angular velocity 
(and the surface density; see below) are non-axisymmetric at leading order in e. 

These expansions may then be substituted into the dynamical equations. First, equation (^l|) at 0(e s ) yields 

- p rQ 2 = -po$', 



which determines the orbital angular velocity. The epicyclic frequency n(r) is defined by 
k 2 = 4Q 2 + 2rf2fi', 

and the dimensionless epicyclic frequency is k — k/Q. The remaining equations may be divided into two sets. 



(66) 
(67) 



4.2 Set A equations 

Equation @ at 0{e s ): 

Equation at 0(e s+2 ): 
(side/, — ~^c) Po = -^-d(Vei. 
Equation @ at 0(e s+1 ): 

po (ttdj, - v n ~ 2p fi(w,£i + rv r u' cos/3) = -{0 cos <f> + 7' sin /? sin </>)<9 c p + (pbo + §po)^<9 c t;ei 

+ — + {ft' cos + 7' sin /3 sin </>) 2 d^(podQV r i) + sin — 7' sin /3 cos <p)d^po- 
Equation @ at 0(e s+1 ): 

po ^09^ — ~~^<J [ v oi + cos</> + 7' sin /3 sin 0)] 

-p fi [r^C + rv rl {0 sin</> - 7' sin /? cos </>)] = -9< po + (pw) + §Mo)^<9 c «ei 

+ + (/3' cos + 7' sin/3sin </>) 2 9^ {poc*c + ru >-i (ft' cos + 7' sin (3 sin 0)] } 
—rQ,(j3' cos + 7' sin /3 sin (/>)(/?' sin </> — 7' sin /3 cos </>)c\po. 
Equation @ at 0(e s+1 ): 



(68) 



(69) 



(70) 



PO (fldj, — ~~^<J (Vfi + J-Vrlf' cos (3) + 
+rfi {ft cos + 7' sin (3 sin 4>)d(po- 



20, 



— + (0 cos + 7' sin/? sin cj>) 2 c\ [po<9c( v 0i + r "rt7' cos/3)] 



(71) 



(72) 



4.3 Set B equations 

Equation © at 0(e s+1 ): 
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-d,/, ) p = ^d^vei - p 



—zd r {r v T i) 8^62 + -d<f,V4,i 



— dr(r 2 Vrl) 8(V02 H — d^v^i 



(0.84, — — d() Pi + [Vrldr ~ — dt + 

V r ) \ r r 

Equation (|3Sj) at 0(e s+3 ): 

(fld<f, — ~dc) pi + (v r idr — -y-d( + ~~®<l>) po = -^-d(v ei — Fpo 
Equation (jZTJ) at 0(e s+2 ): 

PO (fld<p - ~~^J Vl r2 + PO (vrldr - -^c\ + ~~T"^<H V rl + Pi ( ~ Wrl ~~ ~7 t" 91 + rv rl{P' COS + 7' sin /9 

— 2p fi [— \ r ^C + v <j>2 + r(7 + «r27') cos/3 — rv r i{0' sin</> — 7' sin /3 cos </>)(] — ^7(^1 + nvi7' cos /3) 2 
-2pir2(u <# , 1 + rv r ij' cos 0) = -(/3' cos</> + 7' sin /3 sin </>)9 c jpj - (p b o + §£to) 



1 2 1 1 

— 3 r (r « P i) — -d(Vg 2 + -d^v^i 



+{phi + jsVi)-d(Vei\ - {Or - 7' cos 84) po + (pbo + lpo)^dcvei 



+ 



— + (/?' cos0 + 7' sin/3sin</>) 2 c\(poc\u r 2 + pid^v r i) + (/3'cos</> + 7' sin/? sin 0)9^ [^o(9r — 7' cos /3 9^)tvi] 

r 



+—^(d r — 7' cos (3 84) [por 2 (0' cos + 7' sin (3 sin (f)d(V r ij — "" rL (/?' cos + 7' sin /3 sin 4>)8qPo 
+ ~¥ [(^ r - T' cos/3 90)(/i O '" 2 )] c\ [wfli + rv rl (/3' cos</> + 7' sin /3 sin 0)] 



— (dcPo)(dr — 7' cos 84) + v r i(/3' cos + 7' sin /3 sin 0) 

L r 

— (<9cpo)9<£ [(/?' cos + 7' sin ,0 sin </>) (i>0i + ry r i7'cos 0)] + fi'd^po 

+ -(/3' cos + 7' sin /3 sin (j>) (84, I4))dc( v 4>i + rv rij' cos /?) + 0(/3' sin cf> — 7' sin /3 cos <p)8^p,i 

e 8+2\. 



Equation (U2h at 0(e 



po ( — 9 C ) r«6i2 + r(/3 + Wr2/S') cos + r (7 + ^27') sin sin < 



+ P0 (vrldr — "^-^C + ~~^J [ V 81 + rV rl (/#' COS </> + 1 Sm /3 SU1 0)] 

+pi ^9^, — [ u ei +r« r i(/3'cos</> + 7'sin/3sin(/>)] + ^°" rl [v g i + rv r i(0' cos </> + 7' sin /3 sin 0)] 

— pof2 [(w^i + ™ r i7' cos /3)C + ^(/? + v r2 0') sin — r (7 + f r 27') sin cos <f>] 

— — (v^i + rv r ij cos 0) \rQ.C, + rv r i(0' sin 1/1 — 7' sin cos </>)] — [rfi^ + rv r i(0' sin 1/1 — 7' sin /3 cos (/;)] 



-9< |P1 - (PbO + |/io) 



-i-9 r (r 2 w r i) - -d(,ve2 + - 



+ 



+ (0' cos(/> + 7' sin/3sin(/>) 2 9 ? {^o^c [ u »2 + ?"(/3 + v r2 0') cos</) + r(-y + ^27') sin sine/)] 
+ rv r i(0' cos</) + 7' sin sinc/>)] } 
+(/3' cos </> + 7' sin/3sini/>)9c {^(^r — 7' cos /3 9^) [vei + rv r i(0' cos ^ + 7' sin sin 0)] } 
+ -2 (c?r — 7' cos /3 84,) {por 2 (0' cos + 7' sin /3 sin 0)9,; + rv r i(0' cos + 7' sin /3 sine/))] } 
— - (0' cos + 7' sin /3 sin 0) (9cMo) [«ei + »"Wr-i (/3' cos <j> + 7' sin /3 sin </>)] ^ (d(;Vri ) (9 r — 7' cos 84,) (por 2 ) 



+ ^(d( ! po){dr — ^' cos d4,)v r \ + -^{d,;po)d4 > (v4,i + rn r i7' cos/3) 



— -^{84,^0) d((v 4,1 + rvrij' cos/3) — (0' sin</> — 7' sin /3 cos (f>) (0' cos(f> + 7' sin sin </>)9 c [^o("0i + riviV cos/3)] 

— rfi(/3' sin (f> — 7' sin /3 cos </>)(/3' cos </> + 7' sin /3 sin <p)d(pi j(9r — 7' cos /3 9^) [por 3 Sl(0' sin 1/1 — 7' sin /3 cos (/;)] 

—po(0' sin </> — 7' sin /3 cos </>) [(rf2)' + (/3' cos + 7' sin sin </>)c\( u <^i + rn r i7' cos/3)] . 
Equation Q at 0(e s+2 ): 

po ^fi9^, — -^"9c^ \_V42 — |^ri^ 2 + r(-y + u r 27' cos/3) — rv rl (0' sin</> — 7' sin /3 cos </>)(] 

+po ( fri9 r - — Q c 4. ) ( w , a + ru r i7 ; cos/3) + pi ( ri9 - — 9 C ] + r« ri 7' cos/3) + ^j-tV2 + ^^-"ri 

Vrr/ \ r J l\l l\l 
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+—v r i(v ( pi + rv r i"/' cos 0) + — [vei + rv r i (0 cos + 7' sin sin 0)] [rfi£ + rv r i(0' sin 4>~ "/' sin /? cos 0)] 
= —-d$ p + (Hho + §Mo)~dfU0i +rf2'(/3'cos^) + 7'sin/3sin0)a c /ii 

+ ^d r [mo J" 2 (Vfl)'] — rO' (9^po)7' cos/? + (/?' cos</> + 7'sin/3sin</>)c\ [^o(<9 r — j cos0d<f,)(v,i > i + ru r i7' cos /?)] 
+ ^-(9r — 7' cos/3 9^) [/ior 2 (/3' cos</> + 7'sin/3sin</>)c\(««pi + rtvi7' cos/3)] 

+ + (/3' cos</> + 7' sin/3sin</>) 2 9^ {^oc\ [~l r ^C 2 + u <*2 + r(-y + v r 27') cos/3 — rv r i(0' sin — 7' sin /3 cos </>)(] 

+^i<9c(^i +™ r i7'cos/3)} — —d T (fJ,or) ~ -(v</,i + rvrij' cos /?)(/?' cos + 7' sin /? sin 4>)d ( po 

+—CdcHo + i(d<4u r i)(/3'cos</> + 7'sin/3sin</>)(9 c ^o — -(d<hUo)((3' cos + 7' sin /3 sin <f>)dcv r i 
r r r 

a(fyf i °)d<l> + ru r i (/?' cos + 7' sin /? sin 0)] + — (9^o)9< [wei + riVi (/?' cos + 7' sin /3 sin </>)] 

+(/?' sin — 7' sin /3 cos 0) (/3' cos + 7' sin j3 sin </>)c\ {/io [ufli + rv r i (0 cos + 7' sin (3 sin 0)] } 
+A i o(/?' sin </> — 7' sin /3 cos </>)(/?' cos </) + "/' sin /3 sin </>)<9<; [wsi + rv r i{0 cos + 7' sin /3 sin 0)] 
—porD.(0 sin </> — y' sin /3 cos </>) 2 . 



(77) 



4.4 Integrated quantities 

Finally, equation (|^) is also required at 0(e s+2 ), but only in its integrated formjf] 



= 0. (78) 



d T ' yj >J por&4>&Qj + ~d r r j j {pov r2 + piv rl ) r d<j> d( 

The surface density E(r, T) at leading order in e is 

E= J PordC, (79) 

and is independent of </> by virtue of equation (|6^). Other vertically integrated quantities are non-axisymmetric, however, and 
these will be written with tildes. The corresponding azimuthally averaged quantities, written without tildes, are defined by 
the operation 

= (80) 
Thus a suitably averaged radial velocity v r (r,T) may be defined by 

Y,v r = f={F), (81) 
where 

(82) 
(83) 



T — I (pQV r 2 + piVrt) r d( 



is the radial mass flux at 0(e s+3 )J^] Then equation (fF^), divided by 2tv, reads 



a T S + -d r (rEv r ) =0, 



which agrees with equation (g|). Two more definitions will be useful. A suitably averaged kinematic viscosity v(r,T) may be 
defined by 



uT, = V = (V) 
where 



V= / p rd( 

is the vertically integrated viscosity. Finally, the second vertical moment of the density is 



(84) 
(85) 



t Throughout this paper, integrations with respect to tj> are carried out from to 27r, and integrations with respect to £ are carried out 
over the full vertical extent of the disc. 

t There is a radial mass flux at 0(e a+2 ), but its azimuthal average vanishes by virtue of equations (Jt2) and $5^). 
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Z= Por 2 ( 2 rd(, 



(86) 



and its azimuthal average is 



(87) 



4.5 Formal manipulations 

It is clear that the problem under consideration is much more complicated than the analysis of an unwarped viscous disc. 
In the absence of a warp, all quantities are independent of 4> an d symmetric about £ = 0. Furthermore, the quantities v r i, 
V91, vm and vei all vanish. It is worth noting that if the same assumptions could be made for a warped disc, equations (Q) 
and (g) with v\ = v-z = v could be derived directly by integrating equations (fnj) and (0). However, these assumptions 
are inconsistent with the equations of Set A. For example, in equation (^) there is a radial force resulting from the lateral 
imbalance of pressure and viscous stress in a warped disc, and this drives the horizontal velocities v r \ and v$i. 

It is important to note the formal structure of the problem. Set A consists of five coupled non-linear partial differential 
equations (PDEs) with two independent variables {<f>, £} and seven dependent variables {po,po, Po, Pbo, v r i, vei, v<f,i}. The equa- 
tions must be closed by specifying suitable prescriptions for the viscosities and for the thermodynamics. In general, they must 
be solved numerically. Set B is a set of five coupled linear PDEs for the higher-order quantities {pi,pi, pi, Phi, v r2, V92, 1^2}, 
with coefficients that depend on the solutions of Set A and their radial derivatives. Fortunately, all the information required 
from Set B can be extracted by integration, so that only Set A need be solved in detail. 

The aim is to extract equations resembling equations (|5l|)-(l5s|). To this end, multiply equation ^n\j by r, integrate 
vertically and average azimuthally to obtain 

Eu r (r 2 fi)' = h, 

after an integration by parts and the use of equations ( pq ) and (|73|), where 

h ■ 



(^d r [por 4 Q.' — por S + rv rl ^' cos/3) + fi r 3 '(/?' cos^ + 7' sin/? sin (p)d((v^,i + rv r ij cos /?)] 

— po lyei + rv r i{0 cos 4> + 7' sin (3 sin </>)] [rf2£ + rv r i (/?' sin <f> — 7' sin (3 cos 0)1 — po(f3' sin <f> — 7' sin f3 cos 4>)d(;Vri 
+p,or(0' cos <f> + 7' sin f3 sin <f>){f3' sin <j> — 7' sin f3 cos <j>)d( \yei + rv r i(/3' cos <j> + j' sin (3 sin 0)] 

— fior 2 £l((3' sin (f> — 7' sin /3 cos 0) 2 ^ r d£. (89) 

Then multiply equation ( |7rj ) first by (— rsin$>), then by rcos<j!>, in each case integrating vertically and averaging azimuthally 
to obtain first 

Er 2 fi(/3 + v r /3') =h, 

and then 

■Er 2 Sl(-y + v rl ')smf3 = / 3 , 
where 



(90) 
(91) 



I2 = / (r 2 f2^cos0 



1 



d r (r po^ri) + -d^ipov^i) 



H — — d r {r 3 poVri [vei + rv r \{0 cos cf> + J sin ,3 sin 0)] } 

po"c/)i cos [f ei + rivi (/3' cos ^ + 7' sin /3 sin 0)] — 2porQ(^ sin 0(-U0i + rv r ij' cos /3) 
PoTVri sin 4>(v<f,i + rVri'y' cos /3) (/3' sin ^ — 7' sin /3 cos 0) 

sm _ y cog ^g^) {por 2 (/3' cos + 7' sin /3 sin0)9^ [ugi + rtvi(/3' cos + 7' sin /3 sin 0)] } 



-po sin <^(/3' cos + 7' sin (3 sin t^>)c3^ [t>ei + rivi (/?' cos ^ + 7' sin f3 sin < 



+ 



-Or(por d(Vrl) 



fl COS ( 



-d r \por Z 0,{(3' sin — 7' sin /3 cos 0)] + po?" 2 f27' cos /3cos 0(/3' sin <f> — 7' sin /3 ( 



+/io^ sin </>(/?' sin <j>— 7' sin /3 cos </>) [(rfi)' + (/?' cos + 7' sin /3 sin 0)9j(u0i + r« r i7' cos /3)] ^rd( 

and 

I3 = I ^r 2 fi£sin(/> 



(92) 



1 2 1 

-nd r (r pov r x) + -d r f,(pov^ >1 ) 



COS ( 



■9r {r 3 po«ri [wei +r-n r i(,9'cos0 + 7'sin/3sin0)] } 



—poVc/,1 sin [vei + rv r ± ((3' cos (f> + ■y' sin /3 sin </>)] + 2porf2C cos <t>(v<i>i + rv T \^' cos /3) 
+porw r i cos 4>(v,f,i + rv r ij' cos /3)(/3' sin 1^ — 7' sin /3cos <?!)) 
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+ C ° S - (d r — 7' cos/3 9^) {por 2 (/3' cos <j> + 7' sin /3 sin0)9 ? [vgi + rv r i(/3' cos + 7' sin /3 sin0)] } 

+fio cos cos^) + 7' sin /3 sin 0)5^ [v$i + rv r i((3' cos</> + 7' sin /3 sin 0)] — ^^ 9 r (A*or 2 8{i>ri) — ^— — — —d^v^i 

— C ° S - i9 r [/io?" 3 f2(/3' sin </> — 7' sin /3 cos 0)] + po^ 2 !^' cos /? sin 0(/3' sin — 7' sin /3 cos 0) 

— piorcos0(/3'sin0 — 7'sin/3cos0) [(rfi)' + (/?' cos 4> + 7' sin j3 sin c/>)d((v^i + rv r i"f' cos /?)] )rd£. (93) 
Again, it will be useful to consider the complex combination 

h + il3 = J 4r(<9r + 17' cos /?) (V ^ {por 4 £l(,Vri — ipor 3 ivi + rtvi (0 cos + 7' sin /3 sin 0)] 

+ipo7" 3 (/?' cos + 7' sin /3 sin </>)<9c [wei + rtvi {0 cos </> + "/' sin /? sin 0)J — i/xor 2 9^u r i 
— i/j,or 4 Q(P' sin (j> — y' sin /3 cos 0)} ^ + ^e"* {— po(r 2 £l)' (v,-i 

—po(v<t,i + rv rl -y' cos /3) [v B1 + rv r i (/?' cos + 7' sin /? sin <£)] + iporf2C( u </>i + rw r i7' sin /?) 

+ipo»"Wr-i(«d>i + rvrij' cos f3)((3' sin — 7' sin /? cos 0) — —dAvsi + riviV cos/3) 

r 

—iHor(P' sin — 7' sin /3 cos 0) [rfl' + (/?' cos + 7' sin /3 sin (j>)d((v^,i + 7TV17' cos /?)] } ) rd(. (94) 
5 SEPARATION OF VARIABLES 

In order to proceed, it is helpful to make two simplifying, but physically reasonable, assumptions. First, the viscosity coefficients 
are assumed to be locally proportional to the pressure, so that 

p = ap/fl, (95) 

p b = abp/£2, (96) 

where the dimensionless coefficients a(r) and ab(r) are prescribed functions of radius. Secondly, the fluid is assumed to be 
locally polytropic, with T(r) > 1 being a prescribed function of radius. These assumptions allow the equations of Set A to be 
solved by separation of variables. (An isothermal assumption, r = f , would allow a similar simplification.) 

Given that the pressure and density locally satisfy a polytropic relation, p — Kp T , introduce the enthalpy 



dp 



Then equations (p^) and ( |69| ) may be replaced by 

( nd ,- V Jl d< ) ho = ^Z^l d<Vei , (98) 
V r J r 

and Set A is reduced to a problem in four dependent variables {ho, tVi> v$i, {v$\ + rtViV cos /?)}. Note that tp = \ip\ c lx occurs 
only in the combinations 

r(/3' cos <j> + 7' sin/3sin (f>) — \tp\ cos(^> — \), (99) 
r(/3' sin <j) — 7' sin (3 cos <j>) — sin(<^> — x)- (100) 

The solution of Set A is then of the form 

h = r 2 tf [/iW>- X )- |/2(0 -X)C 2 ] , (101) 

t!ri = rn/ 3 (0-xK, (102) 

W9i = rn/ 4 (^-x)C. (103) 

t^i +rw r i7'cos/3 = rn/ 5 (</>-x)C (104) 

where the dimensionless functions fi, - ■ ■ ,fs (whose parametric dependence on r and T has been suppressed) satisfy the 
non- linear ODEs 

ti(4>) = (r-f)/ 4 (0)/iW>), (105) 

ti(<P) = (T+l)/4(0)/a(^), (106) 
/3O/O = hWMt) + *&(<!>)+ [1 + ("b + S")/*^)] /2(0)|^|cos0-a/2(^)/3(0)(l + |^| 2 cos 2 0) 

-a/ 2 (>)|V>|sin0, (107) 
/lW>) = cos<^ + 2/ 3 (^)|V| sin^ + / 4 (0) t/ 4 (0) + cos<^] + 1 - [l + (a b + |a)/ 4 (0)] /2WO 
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-Q/2W [UW + fa(4>)\1>\ cos0] (1 + W 2 cos 2 tj>) + a/ 2 (0)|^ 1 2 cos sin 



+ i(4-K 2 )a/ 2 (^)|V!cos<^, 



(108) 
(109) 



subject to periodic boundary conditions fn{2n) = /n(0). There are five dimensionless parameters S 2 , F, a, «b}- Note that 
equation (105) is decoupled, and that /i admits an arbitrary multiplicative constant which allows the surface density to be 
fitted. The (upper) surface of the disc is given at leading order by 

1/2 



C = C.(0 



/ 2 (4>-x) 



(110) 



The expressions of Section 4.5 may now be simplified. The vertical integrals introduce the quantities I and V which 
are related by 

V = a/ 2 (0- X )fil (111) 

This follows from the definition of V after an integration by parts. Then equations ( p9[ ) and ( |5l"| ) can be identified provided 
that 



?2 = 7772 (/e [(/4 + /3M cos<^)(l + / 3 |<0| sin<£) + af 2 fa\i>\ sincf)- af 2 {fi + f3\ip\ cos ^)|f/'| 2 cos <?!>sin 



+a/ 2 |?/>j 2 sm 2 0] ), 
in which /„ stands for /„ (<fi), and 

contains the azimuthal dependence of Z. To evaluate this, note that 



3/2 



and so 

/«(*) = [/i(*)] 1/(r-1) 



/i(0) 



3/2 



i/(r-i) 



3/2 



Then it follows from equations (105) and (106) that 
f e (4>) = -2/4(0)/ 8 (^), 

while the definition of /e requires (/e) = 1. Note that, for any function F, 

<e*Ffo - X )> = <e iW+x) F(0)) = ^L^F^)). 

m 

Thus 72 + i/3 can be identified with X + iY, provided that 
Qi = T^i^h {-^.h -/*(/» + /sMca»0) +i/ 5 +i/ 3 / 5 |V|sin0- a/2/5 
-ia/ 2 [-i(4- K 2 ) + MV| cos0]|V| sin^} ) 



(112) 

(113) 
(114) 

(115) 

(116) 
(117) 
(118) 



r(e' /6 [/ 3 -i/M/i + /3|^| cos0) + ia/ 2 (/4 + / 3 |V>I cos 0)|V>| cos cj> - ia/2/3 - ia/ 2 |^| sin 0] ) , 



(119) 
(120) 



These correspondences are equivalent: equations (112) and (119), and equations (113) and (120), can be shown to agree 
by using equations (109) and (117). Therefore equations (112) and (120) for Q\ and Q4 are sufficient and will be taken as 
definitive. This completes the demonstration that the fully non-linear problem can be reduced to one-dimensional conservation 
equations for mass and angular momentum as anticipated in Section 3.2. 

In the inviscid case a — Qb = 0, the functions /1, f 2 , fs and f§ are even, while /b and /4 are odd. Then Qi = Q2 = 0, 
i.e. Q4 is purely imaginary. 



6 EVALUATION OF THE COEFFICIENTS 



In principle, equations ( 105 )— ( 109 ) and (117) can be expanded in powers of to determine the dynamics to any desired order. 
In the Appendix, truncated Taylor series for the coefficients Q\, Q 2 and Q3 are developed in this way. This generalizes the 
linear theory of Papaloizou & Pringle (1983) to allow for an arbitrary rotation law, and extends it into the weakly non- linear 



© 1998 RAS, MNRAS 000, 



14 G. I. Ogilvie 




Figure 2. Contour plot of the coefficient Q3 for an inviscid disc with T = 5/3. The horizontal coordinate is the amplitude of the warp. 
The vertical coordinate is the square of the dimensionless epicyclic frequency. The dashed lines bound the region in which a solution is 
found. 



domain. These series can be used to estimate at which amplitude, and in what way, the linear theory breaks down. This 
behaviour is interpreted further in Section 7. The expansion fails only when both k 2 — > 1 and a — > 0; this is the previously 
identified resonant case, which cannot be described usin g this metho d. 

It is also straightforward to integrate equations ( 106 )-( 109 ) and (117) numerically, imposing periodic boundary conditions 
and the normalization condition (fe) = 1. The coefficients Qi, Q2 and Qs can then be determined directly for any amplitude 
of warp, resulting in a fully non-linear theory. Two important cases have been investigated in this way. 



6.1 Inviscid, non-Keplerian disc 

The first case is that of an inviscid disc with F — 5/3 and a — — 0. The remaining parameters are |V>| and k 2 . A contour 
plot of the only non-vanishing coefficient, Q3, is shown in Fig. 2. For reasonably small values of \if)\ there is good agreement 
with the truncated Taylor series. Where k 2 > 1, Qz is negative and the solution can be continued to large values of \ip\, 
although the disc becomes highly non-axisymmetric in this limit. Where k 2 < 1, Q3 is positive and the solution exists only 
for sufficiently small The solution terminates when /2 becomes zero at some point, with the result that the disc can no 
longer contain itself vertically and ruptures. 



6.2 Viscous, Keplerian disc 

The second case is that of a viscous, Keplerian disc with k 2 = 1,F = 5/3 and Qb = 0. The remaining parameters are \%p\ and 
q. Contour plots of the three coefficients are shown in Figs 3, 4 and 5. Again, for reasonably small values of there is good 
agreement with the truncated Taylor series, except when a is small. The solution can be continued to large values of | -0 1 for 
any value of a, and again the disc becomes highly non-axisymmetric in this limit. There are several features to note. First, as 
predicted by the Taylor series, Q\ increases with increasing |i/>| and even becomes positive for sufficiently small a. Secondly, 
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Figure 3. Contour plot of the coefficient Qi for a viscous, Keplerian disc with T = 5/3 and ct\, = 0. The horizontal coordinate is the 
amplitude of the warp. The vertical coordinate is the dimensionless viscosity parameter. 

Q2 has a strong peak at \ip\ — a — where the resonance is located, but the resonance does not extend to large values of 
Thirdly, Q3 is typically much smaller in magnitude than Q2- 



7 SUMMARY AND INTERPRETATION 
7.1 General remarks 

In this paper, the non-linear fluid dynamics of a warped accretion disc has been investigated by developing a theory of fully 
non-linear bending waves for a thin, viscous disc in a spherically symmetric external potential. It has been found that the 
dynamics is described by an equation for the surface density, 

™ + ±JL(rtt, r ) = o, (121) 
at r or 

and an equation for the angular momentum, 

£(Er 2 fi£) + ~^v r r 3 Q£) = ~ (QiTr 2 Q 2 t) +~( Q 2 Xr 3 fi 2 ^) + ±|- ( Q 3 Ir 3 tfe X^)+T. (122) 
at r or r or x ' r or \ or J r or \ Or J 

Here E(r, t) is the surface density, v r (r,t) is the mean radial velocity, fi(r) is the orbital angular velocity, £(r, t) is the tilt 
vector, I(r, t) is the azimuthally averaged second vertical moment of the density [cf. equation (M)] , and a new term T 
represents any additional torque due to self-gravitation, radiation forces, tidal forcing, etc. The dimensionless coefficients Q\, 
Q2 and Q3 depend on the rotation law and the shear viscosity, and also, in the non-linear theory, on the adiabatic exponent, 
the bulk viscosity and the amplitude of the warp. They have been determined both analytically, as truncated Taylor series 
in the amplitude of the warp, and numerically, by solving a set of ordinary differential equations. The angular momentum 
equation can be decomposed into an equation for the component parallel to £, 



E ^^) =11. ( QlJ r 2 tf) - Q 2 lr 2 tf 
dr r or y 
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Figure 4. Contour plot of the coefficient Q2 for a viscous, Keplerian disc with T = 5/3 and «b = 0. 
and an equation for the tilt vector, 

Er 2 nfg+« r ^)=Q 1 Trfi 2 ^ + i|-fQ 2 Tr 3 n 2 ^)+Q 2 Tr 2 n 2 ^ V I£ (Q 3 Zr 3 2 £ x -t x (/ x T). (124) 
V at ar / or r or \ dry or r or \ Or J 

This scheme is a generalization of the form proposed by Pringle (1992), and is equally suitable for numerical implementation. 
The internal torques between one ring in the disc and its neighbours are of three kinds. Coefficient Q\ represents a torque 
tending to spin up (or down) the ring. This would be the usual viscous torque proportional to dtt/dr in a flat disc, but in 
a warped disc there is an additional contribution, not proportional to dfi/dr, due to a correlation between the radial and 
azimuthal velocities induced by the warp; this vanishes in an inviscid disc because the radial and azimuthal velocities are 
perfectly out of phase. Coefficient Q2 represents a torque tending to align the ring with its neighbours, which acts to flatten 
the disc; this also vanishes in an inviscid disc.|] Coefficient Q3 represents a torque tending to make the ring precess if it is 
misaligned with its neighbours; this leads to the dispersive wave-like propagation of the warp. It is likely that this form of 
angular momentum equation applies in much more general circumstances than those considered in this paper. However, the 
present analysis provides a consistent evaluation of the three coefficients under a number of simplifying assumptions. 

There are essentially four different assumptions or approximations involved in this analysis. The first assumption is that 
the disc is thin (H/r <C 1), which is essential for the asymptotic expansions in Section 4. The second assumption is that 
the fluid obeys the compressible Navier-Stokes equation even though, in reality, the turbulent stresses in an accretion disc 
are likely to be anisotropic and may not be purely viscous. The third assumption is that the fluid is locally polytropic and 
has viscosity coefficients locally proportional to the pressure. These are the simplest possible closure relations in the absence 
of a full treatment of thermal and radiative physics, and greatly simplify the problem by allowing separation of variables 
(Section 5). In a fully self-consistent treatment, the thickness of the disc, and therefore the relation between X and E, would 
be determined by the local viscous dissipation of energy. The final assumption is that the conditions for resonance, 

§ It is to be expected that Q2 will normally be positive, since it represents a diffusion coefficient. The linear theory (see the Appendix) 
suggests that this is true except possibly for unphysical rotation laws (k 2 > 8). 
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Figure 5. Contour plot of the coefficient Qz for a viscous, Keplerian disc with Y = 5/3 and «b = 0. 

" 2 <,H/r and a<,H/r (125) 

are not simultaneously satisfied. (It has been shown in Section 6.2 that a further condition for resonance is that the amplitude 
of the warp be small.) The asymptotic expansions formally break down as the resonance is approached. The nature of the 
approximation is as follows. In equations (|7^) and (f?^) governing the horizontal velocities induced by the warp, the time- 
derivatives of the velocities (in the inertial frame) do not appear because those terms are of higher order in the assumed ordering 
scheme. The horizontal velocities are in a 'geostrophic' balance in which the driving forces are matched instantaneously by 
inertial and viscous forces, and for this reason they do not constitute additional dynamical degrees of freedom. In the resonant 
case, this approximation breaks down and the time-derivatives must be restored. This results in a different ordering scheme 
and makes the problem hyperbolic rather than parabolic (Papaloizou & Lin 1994). The non-linear dynamics of the resonant 
case is likely to be very different. 

Two important cases merit further interpretation. 




7.2 Inviscid, non-Keplerian disc 

In the inviscid case, the radial velocity vanishes, the surface density is independent of time, and the dynamics is described by 
an equation for the tilt vector of the form 

where Qs is plotted in Fig. 2 for the case V = 5/3. It has the Taylor series (see the Appendix) 

Qs = + 4( 3_yi ia)3 l^ + (127) 

where 
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= r 



dl 

dr 



is the amplitude of the warp. Equivalently, 

lr 3 Sl 4 (6 + Y)lr°Sf 



at 



dr 



2(Sl 2 



) + 4(3-r)(fi 2 - k 2 ) 2 



+ o 



oe 4 

dr 



I X 



01 

dr 



(128) 



(129) 



To make the connection with the linear theory of bending waves, suppose that the disc is initially flat, with L — e z , and is 
then perturbed weakly, with \£ x \, \£ y \ <C 1. Then 



v 2 o dW 



ld_ 

r dr 



2r 3 n« 



2(Sl 2 ~ k 2 ) 



or or or 



+ 



(6 + r)Jr 5 n 6 

4(3-r)(n 2 -k 2 ) 



;i\Wr^+0(W S ) 



(130) 



where W = £ x +i£ y . This is a Schrodinger-type wave equation with non-linear dispersion. The linearized form of this equation, 



ld_ 

r dr 



Xr 6 Sl 



3n4 



f.dW\ 



2(Sl 2 - k 2 ) 



(131) 



may be compared with equation (12) of Papaloizou & Lin (1994). Taking the low- frequency, non-resonant limit of their 
equation, setting the precession frequency to zero for a spherical potential, restoring the time-derivative and identifying 
W <-+ g* , one obtains 



at 



dr 



1SI' 



2(0? - 



k 2 ) V dr ) 



(132) 



The only discrepancy concerns whether a factor r SI should appear inside or outside the radial derivative. This is of little 
significance for nearly Keplerian discs in which r 3 Sl 2 is approximately constant, and may be attributable to a failure of the 



approximations adopted in the passage from equation (11) to equation (12) of that paper. It is evident that equation (131) is 
of the correct form to ensure the conservation of angular momentum. 



In a short-wavelength (WKB) limit, equation (13C) has wave-like solutions of the form 



W(r,t) = W(r)exp 



icot + i / k(r)Ar 



(133) 



where u is the frequency, k(r) is the radial wavenumber (satisfying \kr\ 2> 1) and W(r) is a slowly varying function. The 
non-linear WKB dispersion relation is then 



- = ± 1 > 
SI 2 \ SI 2 



Ik 2 



i + 



(6 + r) 
2(3 -r) V SI 2 



n 2 



2 \w\ 2 + o(\w\ 4 



(134) 



the ± arising since the complex-conjugate solution is equally valid but is physically distinct. In the linear approximation, this 
can be shown to agree with the analytic WKB dispersion relation for isothermal discs found by Lubow & Pringle (1993). 
To make the connection, note that their K x is our kH , where H is the isothermal scale height of the disc, while their F 
is our (ui ± SI) /SI for an m = =pl mode. Expanding their equation (54) about the point K x = 0, F = ±1, and assuming a 
non-Keplerian disc, one obtains for the 'n = 0' mode 



uj i f sr 

Sl~ 2 I SI 2 — k 2 



(kH) 2 + O ((kH) 4 ) 



(135) 



For an isothermal disc X = T.H 2 and the two agree. 

The interpretation of the non-linear term in the dispersion relation is that the dispersion of waves either increases or 



decreases as the amplitude increases, depending on the signs of (SI 



and of (3 — r). The cubic non-linearity arises 



through a three-mode coupling involving (i) the 'tilt' mode or f mode (m = 1, of odd symmetry), consisting locally of a 
uniform vertical translation of the disc, (ii) an inertial or r mode (m = 1, of odd symmetry), consisting of a horizontal 
epicyclic motion proportional to £, and (iii) an acoustic or p mode (m = 2, of even symmetry), consisting of a vertical motion 
proportional to f. The magnitude of the coupling depends on the adiabatic exponent because of the compressive nature of 
the p mode, which has a natural frequency of (1 + T) lll2 Sl in the comoving frame. Indeed, in the unlikely case V = 3, the p 
mode would be driven at resonance and a large response would result. 

Fig. 2 shows that, when F = 5/3, the non-linear behaviour is very different depending on the sign of (S1 2 ~k 2 ). If n 2 > SI 2 
the disc can support warps of large amplitude and the dispersion coefficient becomes smaller in magnitude as the amplitude 
increases. However, if k 2 < SI 2 , the dispersion coefficient increases with increasing amplitude and the disc eventually ruptures. 
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7.3 Viscous, Keplerian disc 

In the viscous, Keplerian case, the three coefficients are plotted in Figs 3, 4 and 5 (with V — 5/3 and a b =0). They have 
Taylor series (see the Appendix) 



, 3a 

?! = — + 



1 - 17a 2 + 21a 4 
4a(4 + a 2 ) 



l^| 2 +0(l^i 4 ), (136) 



Q 2 +iQ 3 = 1 9 + ^ + 6 f + { r - + ^ + g [r 2 i(a b + |a)] \^ + om% m 

2a(2 + ia) 1 4a [(3 - F) + 2i(a b + |a)] (2i - a) 3 (2i + a) J 1 <m " V ; 

where 

a = 12 + 115ia + 49a 2 + 109ia 3 - 408a 4 - 41ia 5 + 2a 6 +4ia 7 , (138) 

b = 18 + 87ia-87a 2 -24a 4 -6ia 5 , (139) 

£ = 12 + 25ia - 36a 2 - Ilia 3 + 140a 4 + 21ia 5 + 2a 6 . (140) 

In the important limit a < 1, db C 1 (but requiring a J> H/r to avoid the resonant case), 

Qi « ~ + ^iVf + OWl*), (141) 
Q 2 « ^ + 0(|^| 2 ), (142) 

Q 3 « | + 0(|^| 2 ). (143) 

This correctly predicts the behaviour seen in Fig. 3, in which the usual viscous torque parallel to t can be reversed by a warp 
of amplitude |t/>| Jj \/24a. This occurs because the usual viscous torque is small (proportional to a) and is easily overwhelmed 
by the torque resulting from the correlation of the radial and azimuthal velocities, since these are almost resonantly driven. It 
also shows that the diffusion of the warp is more important than any dispersive wave propagation in this limit. Unfortunately, 
the truncated Taylor series for Q2 is inaccurate even for reasonably small values of in this limit, and does not predict that 
it ultimately decreases with increasing \if)\, as seen in Fig. 4. Therefore a numerical solution is required in this case. 

To the extent that the dispersion coefficient Q3 can be neglected, the original equations of Pringle (1992) are formally 
valid, but with the following caveats concerning the viscosities v\ and V2- For small- amplitude warps, the approximations 

v x « v (144) 

and 



V 2 



2(1 + 7a 2 



a 2 (4 + a 2 ) 



for a < 1 (145) 

2a z 



hold (cf. Papaloizou & Pringle 1983). These are the approximations relevant for deciding whether a flat disc is unstable to 
the radiation-driven instability (Pringle 1996). The fact that V2 is potentially much larger than v\ means that the radius 
outside which the instability sets in is likely to be much larger than would be estimated on the assumption that v% = v\. If 
the instability does proceed, then v\ and 1*2 are subject to non-linear corrections depending on the amplitude of the warp. 
In the case investigated in this paper, both v\ and V2 typically decrease with increasing amplitude, with v\ even becoming 
negative if a is sufficiently small. Although this means that the usual accretion torque is reversed, the disc so strongly resists 
being warped in this limit that this situation may not persist for long unless the warp is strongly forced. 



7.4 Outlook 

This work has shown that, with the exception of discs that are both accurately Keplerian and almost inviscid, the dynamics 
of a warped accretion disc can be reduced to simple one-dimensional conservation equations for mass and angular momentum 
even if the warp is non-linear. This generally supports the approach adopted by Pringle (1992), although, for completeness, 
it requires an additional type of internal torque to be included. Based on the assumption of an isotropic underlying viscous 
process, it also predicts the values of the coefficients in the equations and their variation with local parameters of the disc 
and with the amplitude of the warp. 

Some features of the behaviour of the system can be estimated by inspection of the equations and the variation of the 
coefficients. Non-linear effects are strongest in the case of most interest, that of a Keplerian disc with H/r ^a<l. However, 
the detailed application of this work to astrophysical situations properly requires a numerical solution of the equations. As 
shown by Pringle (1992 et seq.), the reduction of the problem to one- dimensional equations which do not require the fast 
orbital time-scale to be followed offers very significant practical advantages over three-dimensional numerical simulations. 

At the same time, there are two important ways in which this analysis should be improved in future. First, a better 
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analytical modelling of the turbulent stress tensor in accretion discs is desirable. It should be possible, using local numerical 
simulations, to measure the response of magnetohydrodynamic turbulence to imposed motions such as those induced by a 
warp (Torkelsson et al., in preparation) and to use this to calibrate an analytical model. Secondly, the effect of increased 
dissipation caused by the warp on the thickness of the disc, and therefore on the torques between neighbouring rings, ought 
to be taken into account. Unfortunately, to do this properly requires the solution of an energy equation including radiative 
transport, which is unlikely to be amenable to the technique of separation of variables used in this paper. However, rather 
than invalidating the form of the equations derived here, these developments would be expected only to provided improved 
values for the coefficients. 
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APPENDIX A: TRUNCATED NON-LINEAR EQUATIONS 

The aim of this section is to derive truncated Taylor series for the coefficients Qi and Q<t = Qi + 1Q3, which will provide an 
equation for the tilt vector correct to cubic order. To achieve this, introduce the expansions 

M4>) = /10 + W 2 /i2(<rt + o(|^| 4 ), (Al) 

h{4>) = /20 + W 2 /22(<rt + o(|^| 4 ), (A2) 

M4>) = Whi(4>) + W\ka(cb)+0(W 5 ), (A3) 

U(4>) = W 2 f i2 W + o(\i>\ 4 ), (A4) 

H4>) = Wf 51 (4>) + M 3 / 5 3W>) + o(H 5 ), (A5) 

fe(4>) = f 60 + M 2 f 6 2W + O(W 4 ). (A6) 

Note that each function is either even or odd in \ip\, and that all terms with the scaling l^l are axisymmetric, since they 
represent an unwarped disc. 

Al Zeroth-order solution 



At zeroth order, the solution is that of an unwarped disc. For the vertical equilibrium, equation ( 106 ) at O(l) yields simply 
/20 = 1. (A7) 
Also 

feo = 1, (A8) 
since (fe) = 1 by definition. 
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A2 First-order solution 



The horizontal velocities at first order are determined by equation (107) at 0(|i/?|), 
/31 (4>) + a/31 (4>) ~ 2/51 (0) =cos4>-a sin 0, 



and equation (109) at O(|-0| ) , 

/^)+a/ 51 (0) + i£ 2 /3iW>) = i(4-K 2 )acos0. 

The solution is of the form 

/31 (0) = CVi cos + Sri sin 0, 
fel{<t>) = C^i cos0 + S$i sin0, 

with 



-2 
1 

a 

This may be expressed more compactly using a complex notation, with (generically) 
Z = C* + iS, 
so that 



a 


1 


-2 


-1 


a 





2 K 





a 





2 K 


-1 





" CVi " 




1 




Sri 




—0 








i(4-£ 2 )a 




- S,/,! _ 








a-i -2 








1 — ia 


^k 2 a — i 








_i(4-S 2 )a_ 



The determinant of this matrix is 
-(l-K 2 )-2ia + a 2 , 



(A9) 

(A10) 

(All) 
(A12) 

(A13) 

(A14) 
(A15) 
(A16) 



and so a solution exists unless the disc is both Keplerian (k 2 — 1) and inviscid (a = 0) to the accuracy of these scalings. 
(As anticipated, the resonant case cannot be described using the present expansions.) The solution follows by inversion of the 
matrix. In detail, 



Z r l 
Z^l 



i — (4 — k 2 )a + ia 2 
(l-S 2 )+2ia-a i 

1 

2 



k 2 + 2i(2 - k 2 )a - (4 - k 2 )a 2 



(1 - k 2 ) + 2ia - a 2 



(A17) 
(A18) 



A3 Second-order solution 



The enthalpy and vertical velocity at second order are determined by equation (105) at 0(1^1 

/r2«0 = (r-l)/42«>)/lO, 



equation ( |106|) at 0(|i/>| 2 ), 

/22«>) = (r+l)/42«>), 



and equation (108) at 0(1^1 )i 

/42(<W + (Ob + |q)/42(0) = -/3i(0) cos0 + (2 sin <f> - a cos <j>)f 31 (0) - /22(0) +acos<? 
These may be combined to give 

fi2{<t>) + («b + H/«(0) + (r + l)/4 2 (0) = (3CVl - «S r l + Ot) COS 20 + (3S r l + oCVl) 

The solution is of the form 

/42 (0) = Ce2 cos 20 + So2 sin 20, 

with 



sin 20. 



-(3-T) 2(a b + |a) 
-2(a b + |a) -(3-T) 



Cff2 
Sb2 



3C r i — aSri + a 
3<SVi + etCVi 



(A19) 
(A20) 
(A21) 
(A22) 
(A23) 

(A24) 
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In a complex notation, the solution is 
„ —(3 + \d)Z r \ — a 



(3-T) + 2i(a b + | a) 

-3i + 2(6 - ft 2 )a - (1 + ft 2 )ia 2 + 2a 3 
[(3 - T) + 2i(a b + §«)] [(1 - ft 2 ) + 2ia - a 2 ] ' 

It then follows that 

M0) = -|(r- l)/io(5 , S2 cos20-C e2 sin2< ? !>), 

/22WO = -i(r + l)(5e 2 cos20-C , e2 sin20) + i5 r i-iaC , ri, 

/e2 {<t>) = Se 2 cos 2^ - C02 sin 20. 



(A25) 



(A26) 
(A27) 
(A28) 



A4 Third-order solution 



The horizontal velocities at third order are determined by equation (107) at 0(|t/>| 3 ), 

fM + af 33 (4>) - 2/ 53 (0) = / 42 O/O/31 (<£) + [/ 22 (0) + (a* + |a)/ta(0)] cos0 - a [/ 22 (0) + cos 2 0] / 31 (0) 
-a sin 0/ 22 ((/>), 

and equation (]T09|) at 0(|^| 3 ), 



(A29) 

(A30) 

(A31) 
(A32) 

(A33) 

u — 1 J I ^03 J I "2 

where 

ifc = [| - |(r + l)ia] Z e2 ^rt + [|(r + l)(i - a) + |(a b + §a)] Z 92 - |aZ P i - §(&•! - aC r i) [aZri - (1 - ia)] 

+±iaS rt , (A34) 
i? 2 = [5 - i(r + l)ia] Zea-Zji + |(4- ft 2 )(r + l)iaZ fl2 - faZ^i - §a(SVi - aC r i) [Z^i - §(4- ft 2 )] + |ia5^i. (A35) 
The m = 3 terms will not be required. 



fsaW + a/saW + \*hz{4>) = /42(0)/ 5 i(0) - a + cos 2 4>] fia(4) + |(4- ft 2 )acos0/ 22 (0). 

The solution is of the form 

/33 (</>) = C,-3 cos + SV3 sin <f> + {m = 3 terms}, 
/5s(0) = C^3Cos^ + S , 3 sin0 + {m = 3 terms}, 
with 



"a-i -2 




ZrS 




Ri~ 


_ ift 2 a-i_ 




Z<f,3 




R2 



A5 Evaluation of the coefficients 

The coefficients Q\ and Q4 are even in \ip\ 

Qi = Qio + |V>! 2 Qi2+0(|V>| 4 ), 
Qa = Q i0 + W 2 Q42+O(W 4 ). 

At leading order, 

Q10 = -±(4-ft 2 )a 
and 

Q40 = (e"* [(1 - iaO/31 (</>) - ia sin 0] ) 
= ia+Kl-ia)^! 

"i- 2a + (7- ft 2 )ia 2 ~ 



and have expansions 



(A36) 
(A37) 

(A38) 



(1 - ft 2 ) + 2ia - a 2 
At second order, 

Q12 = -|(4 - S 2 )a(5 ri - aC rt ) - i(C ri C 01 + 5 ri S # ) + |aC i 



(A39) 



-(8 - 12ft 2 + 3ft 4 )a - (2 - ft 2 )(28 - 12ft 2 + ft 4 )a 3 + (8 - ft 2 )(4 - ft 2 )a 5 
(l-ft 2 ) 2 + 2a 2 (l + ft 2 )+a 4 



(A40) 



© 1998 RAS, MNRAS 000, 



The non-linear fluid dynamics of a warped accretion disc 23 



and 

Qa2 = (e"* {(1 - ia)f S3 ((/>) - i [/31OA) - acos<j>] [/42(» + /3i(^)cos0] - iaf 22 (</>), f 31 (<t>) - 10/22(0) sin0 
+/ 6 2(0) [(1 - i«)/3i(0) - iasin^]} ) 
= |(1 - ia)Z r3 + i [-i+ |(r - 1)q] + |(3 - r)iaZ 92 - \iZ rl Z* rl - §LZ* a + \\aZ rl 

-\la{S r i - aCri)(Z rl + i) + i Q S r i 

a + br + c[r-2i(g b + |a)] 

4 [(3 - T) + 2i(a b + fa)] [(1 - ft 2 ) + 2ia - a 2 ] 3 [(1 - ft 2 ) - 2ia - a 2 ] ' 

where 

a = 6i(l - ft 2 ) 2 - 3(1 - ft 2 )(8- 7ft 2 + ft 4 )a + 3(l - ft 2 )(28- 32k 2 + 2k 4 + ft 6 )ia 2 

-(258 - 337ft 2 + 103ft 4 - 15ft 6 + 3ft 8 )a 3 + (224 - 750ft 2 + 511ft 4 - 106ft 6 + 6ft 8 )ia 4 

-(442 - 511ft 2 + 116ft 4 + 2ft 6 )a 5 - (552 - 560ft 2 + 124ft 4 - 7ft 6 )ia 6 + (562 - 153ft 2 - ft 4 )a 7 

+ (66 - 26ft 2 + ft 4 )ia 8 + 2(1 - 2ft 2 )a 9 - 4ia 10 , (A42) 
b = -(4- ft 2 )a [(1 - ft 2 ) + 2ia - a 2 ] [(1 - ft 2 ) - (2 - ft 2 )ia + a 2 ] [3 + 2(6 - ft 2 )ia + (1 + ft 2 )a 2 + 2ia 3 ] , (A43) 
c = i(l - ft 2 ) 2 - 2(2 - ft 2 )(l - ft 2 )(l + ft 2 )a - (2 - ft 2 )(l - ft 2 )(2 + 9ft 2 - 3ft 4 )ia 2 - 2(25 - 40ft 2 + 24ft 4 - 3ft 6 )a 3 

-(262 - 428ft 2 + 233ft 4 - 44ft 6 + 2ft 8 )ia 4 + 2(115 - 119ft 2 + 23ft 4 - ft 6 )a 5 

+ (164 - 196ft 2 + 45ft 4 - 2ft 6 )ia 6 - 2(87 - 18ft 2 + ft 4 )a 7 - 3(9 - 2ft 2 )ia 8 - 2a 9 . (A44) 
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